Figure 1 - This figures shows…
library(png)
library(maptools)
Checking rgeos availability: TRUE
library(raster)
library(gam)
Loading required package: splines
Loading required package: foreach
foreach: simple, scalable parallel programming from Revolution Analytics
Use Revolution R for scalability, fault tolerance and more.
http://www.revolutionanalytics.com
Loaded gam 1.14
setwd("~/Desktop/Botero postdoc 2016/Human density and the origins of agriculture/")
par(mar=c(0,0,0,20))
d <- readPNG("Larson_dates.png")
plot(seq(0,18, length.out = 19), seq(0,36, length.out = 19), type="n",ylim=c(0,36),xlim=c(0, 18), xaxt="n")
rasterImage(d, 0,0,18,36, interpolate=TRUE, col=d)
Start_of_early_window <- 16-12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 17-4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
These dates are provided in the supplimentary information for the Larson (2014) paper. I’ve copied those values into a .csv table provided here.
domestication_times <- read.csv("Domestication timing larson 2014.csv")
dim(domestication_times)
[1] 77 8
| Region | Species | Start.Exploitation | Finish.Exploitation | Start.predomestication | Finish.predomestication | Start.Domestication | Finish.Domestication |
|---|---|---|---|---|---|---|---|
| Southwest asia | Wheat | 12.00 | 11.25 | 11.25 | 11.00 | 11.00 | 9.00 |
| Southwest asia | Barley | 12.00 | 11.25 | 11.25 | 10.50 | 10.50 | 9.00 |
| Southwest asia | Lentil | 12.00 | 11.00 | 11.00 | 10.50 | 10.50 | 9.00 |
| Southwest asia | Pea | 11.50 | 11.00 | 11.00 | 10.00 | 10.00 | 8.50 |
| Southwest asia | Chickpea | 11.00 | 10.50 | 10.50 | 10.25 | 10.25 | 8.25 |
| Southwest asia | Broadbean | NA | NA | NA | NA | 10.50 | NA |
| Southwest asia | Flax | 12.00 | 9.50 | NA | NA | 9.50 | NA |
| Southwest asia | Olive | 10.00 | 6.00 | NA | NA | 6.00 | NA |
| Southwest asia | Sheep | 12.00 | 10.50 | 10.50 | 9.75 | 9.75 | 8.00 |
| Southwest asia | Goat | 12.00 | 10.50 | 10.50 | 9.75 | 9.75 | 8.00 |
| Southwest asia | Pig | 12.00 | 11.50 | 11.50 | 9.75 | 10.25 | 9.00 |
| Southwest asia | Cattle, taurine | 11.50 | 10.50 | 10.50 | 10.25 | 10.25 | 8.00 |
| Southwest asia | Cat | NA | NA | 10.50 | 4.00 | 4.00 | NA |
| South Asia | Tree cotton | 8.50 | 4.50 | NA | NA | 4.50 | NA |
| South Asia | Rice | 8.00 | 5.00 | 5.00 | 4.00 | 4.00 | 2.50 |
| South Asia | Little millet | NA | NA | NA | NA | 4.50 | NA |
| South Asia | Browntop millet | NA | NA | NA | NA | 4.00 | NA |
| South Asia | Mungbean | NA | NA | 4.50 | 3.50 | 3.50 | 3.00 |
| South Asia | Pigeonpea | NA | NA | NA | NA | 3.50 | NA |
| South Asia | Zebu cattle | 9.00 | 8.00 | NA | NA | 8.00 | 6.50 |
| South Asia | Water buffalo | 6.00 | 4.50 | NA | NA | 4.50 | NA |
| East Asia | Broomcorn Millet | 10.00 | 8.00 | NA | NA | 8.00 | NA |
| East Asia | Foxtail millet | 11.50 | 7.50 | NA | NA | 7.50 | NA |
| East Asia | Rice | 10.00 | 8.00 | 8.00 | 7.50 | 7.50 | 5.00 |
| East Asia | Soybean | 8.50 | 5.50 | NA | NA | 5.50 | 4.00 |
| East Asia | Ramie | NA | NA | NA | NA | 5.25 | NA |
| East Asia | Melon | 7.00 | 4.00 | NA | NA | 4.00 | 3.75 |
| East Asia | Pig | 12.00 | 8.50 | NA | NA | 8.50 | 6.00 |
| East Asia | Silkworm | 7.00 | 5.25 | NA | NA | 5.25 | NA |
| East Asia | Yak | NA | NA | NA | NA | 4.25 | NA |
| East Asia | Horse | 7.50 | 6.75 | 6.75 | 5.50 | 5.50 | 4.00 |
| East Asia | Bactrian Camel | NA | NA | NA | NA | 4.50 | NA |
| East Asia | Duck | 2.50 | 1.00 | NA | NA | 1.00 | NA |
| East Asia | Chicken | 6.00 | 4.00 | NA | NA | 4.00 | NA |
| New Guinea | Banana | 10.00 | 7.00 | 7.00 | 4.00 | 4.00 | NA |
| New Guinea | Taro | 10.00 | 7.00 | 7.00 | 4.00 | NA | NA |
| New Guinea | Yam | 10.00 | 7.00 | 7.00 | 4.00 | NA | NA |
| Africa and Arabia | Date palm | 7.00 | 6.00 | NA | NA | 5.00 | NA |
| Africa and Arabia | Sorghum | 8.00 | 4.00 | NA | NA | 4.00 | NA |
| Africa and Arabia | Pearl millet | NA | NA | NA | NA | 4.50 | 3.50 |
| Africa and Arabia | Fonio | NA | NA | NA | NA | 2.50 | NA |
| Africa and Arabia | Cowpea | NA | NA | NA | NA | 3.75 | NA |
| Africa and Arabia | Hyacinth bean | NA | NA | NA | NA | 3.75 | NA |
| Africa and Arabia | Rice | 3.50 | 2.00 | NA | NA | 2.00 | NA |
| Africa and Arabia | Oil palm | 9.25 | 3.50 | NA | NA | 3.50 | NA |
| Africa and Arabia | Cattle | NA | NA | 9.00 | 7.75 | 7.75 | 6.50 |
| Africa and Arabia | Donkey | 9.00 | 5.50 | NA | NA | 5.50 | 3.50 |
| Africa and Arabia | Dromedary camel | 6.50 | 3.00 | NA | NA | 3.00 | NA |
| Africa and Arabia | Guinea fowl | NA | NA | 2.50 | 1.50 | 1.50 | NA |
| North America | Squash | 6.50 | 5.00 | NA | NA | 5.00 | NA |
| North America | Sunflower | 6.00 | 4.75 | NA | NA | 4.00 | NA |
| North America | Sumpweed | 6.00 | 4.50 | NA | NA | 4.00 | NA |
| North America | Pitseed goosefoot | 4.75 | 3.75 | NA | NA | 3.75 | NA |
| Meso-america | Squash (pepo) | NA | NA | NA | NA | 10.00 | 9.50 |
| Meso-america | Maize | 10.00 | 9.00 | NA | NA | 9.00 | NA |
| Meso-america | Foxtail millet-grass | NA | NA | NA | NA | 6.00 | 4.00 |
| Meso-america | Common bean | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Avocado | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Chile pepper | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Turkey | NA | NA | NA | NA | 2.00 | NA |
| South America | Chili pepper | NA | NA | NA | NA | 6.00 | NA |
| South America | Peanut | NA | NA | NA | NA | 5.00 | NA |
| South America | Cotton | NA | NA | NA | NA | 6.00 | NA |
| South America | Coca | NA | NA | NA | NA | 8.00 | NA |
| South America | Now-minor root crops (arrowroot, leren) | NA | NA | NA | NA | 9.00 | NA |
| South America | Squash | NA | NA | NA | NA | 10.00 | NA |
| South America | Common bean | NA | NA | NA | NA | 5.00 | NA |
| South America | Lima bean | NA | NA | 8.25 | NA | 6.00 | NA |
| South America | Monioc | NA | NA | NA | NA | 7.00 | NA |
| South America | Sweet potato | NA | NA | NA | NA | 5.00 | NA |
| South America | White potato | 7.00 | 4.50 | NA | NA | 4.00 | NA |
| South America | Quinoa | 5.00 | NA | NA | NA | 3.50 | NA |
| South America | Yam | NA | NA | NA | NA | 5.50 | NA |
| South America | Llama | 10.00 | 6.00 | NA | NA | 6.00 | NA |
| South America | Alpaca | 10.00 | 5.00 | NA | NA | 5.00 | NA |
| South America | Guinea pig | NA | NA | NA | NA | 5.00 | NA |
| South America | Muscovy Duck | NA | NA | NA | NA | 4.00 | NA |
par(mar=c(5,4,6,1))
dates <- unlist(domestication_times[3:8])
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset" )
mtext("This tells us about how evenly our evidence is distributed in time", 3, line=1)
hist(dates, breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset with Larson(2014) date windows")
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)
mtext("Early Holocene", 3, line = -1, adj=.3)
mtext("Middle Holocene", 3, line= -1, adj=.6)
par(mfrow=c(2,3), mar=c(4,4,2,0))
specific_dates <- domestication_times[3:9]
for(i in c(1, 3, 5, 2, 4, 6)){
hist(specific_dates[,i], breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main= names(specific_dates)[i])
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(specific_dates[,i], breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)
}
I’m creating new rows for this table, combining dates in different ways to make the CDFs below look more authentic. This makes it so that pre-ag always happens before post-ag. What I’ve done is given the later date to the earlier date when those dates are missing.
h <- which(is.na(domestication_times[,3]))
domestication_times <- cbind(domestication_times, rep(NA, length(domestication_times[,1])))
domestication_times[,9] <- domestication_times[,3]
domestication_times[h,9] <- domestication_times[h,7]
colnames(domestication_times)[9] <- "adopt exploitation date"
domestication_times[,10] <- domestication_times[,7]
domestication_times[which(is.na(domestication_times[,10])),10] <- 0
colnames(domestication_times)[10] <- "start of ag"
#save(domestication_times, file="~/Desktop/Human density and the origins of agriculture/Domestication timing larson 2014.Rdata")
I think these are best described by a cummulative distribution, showing how they accumulate over time.
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(levels(domestication_times$Region)[ type_number])
print(match)
print(j)
}
[1] "Africa and Arabia"
[1] 7.00 8.00 4.50 2.50 3.75 3.75 3.50 9.25 7.75 9.00 6.50 1.50
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
[1] "East Asia"
[1] 10.00 11.50 10.00 8.50 5.25 7.00 12.00 7.00 4.25 7.50 4.50 2.50 6.00
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
[1] "Meso-america"
[1] 10 10 6 3 3 3 2
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 4, 7, 8
[1] "New Guinea"
[1] 10 10 10
Empirical CDF
Call: ecdf(maxer - match)
x[1:1] = 0
[1] "North America"
[1] 6.50 6.00 6.00 4.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 0.5, 1.75
[1] "South America"
[1] 6.0 5.0 6.0 8.0 9.0 10.0 5.0 6.0 7.0 5.0 7.0 5.0 5.5 10.0 10.0 5.0 4.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
[1] "South Asia"
[1] 8.5 8.0 4.5 4.0 3.5 3.5 9.0 6.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
[1] "Southwest asia"
[1] 12.0 12.0 12.0 11.5 11.0 10.5 12.0 10.0 12.0 12.0 12.0 11.5 4.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:6] = 0, 0.5, 1, ..., 2, 8
par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
#print(j)
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq, ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
if(i == 2 | i == 1)axis(2)
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
}
par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
#print(j)
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq, ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
abline(v= maxer - quantile(j)[2], col="limegreen", lwd=2)
if(i == 2 | i == 1)axis(2)
if(i == 2)mtext("25%", 3, line=3.5, adj=-1, col="limegreen")
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
if(i == 4)mtext("Choose a y to predict an x", 3, line=3.3, col="cornflowerblue")
break_one <- maxer
break_two <- maxer - quantile(j)[2]
polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.2), border=adjustcolor("cornflowerblue",alpha=1))
lines(x=c(break_two, break_two), y=c(0,-1), col="cornflowerblue")
abline(h = 25, col="limegreen", lwd=2)
}
Make this a function. There is a choice of two methods here. At the end of this section we need to print the desision we’re passing to the later analyses.
origins
class : SpatialPolygonsDataFrame
features : 20
extent : -104.7263, 144.9554, -26.31802, 43.00233 (xmin, xmax, ymin, ymax)
coord. ref. : NA
variables : 8
no non-missing arguments, returning NAno non-missing arguments, returning NAno non-missing arguments, returning NAno non-missing arguments, returning NAno non-missing arguments, returning NAno non-missing arguments, returning NA
names : CONTINENT, SQMI, SQKM, OID, CONTINEN_2, SQMI_2, SQKM_2, OID_2
min values : C/S_Andes, NA, NA, NA, Europe, 3821854, 9898597, 7
max values : West Africa T, NA, NA, NA, Europe, 3821854, 9898597, 7
library(maps)
map()
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))
database does not (uniquely) contain the field 'name'.
map()
d <- readPNG("Larson_origins.png")
rasterImage(d, -180, -90, 180, 110, interpolate=TRUE, col=d)
map(add=TRUE)
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))
database does not (uniquely) contain the field 'name'.
This is obviously a bad projection fit right now.
#subset and reorder origins. This is currently done at the end of the plot but should be moved forward.
# Load data for population density
load("PopD_all_December.rdata")
PopD.ALL
class : RasterStack
dimensions : 288, 720, 207360, 18 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -60, 84 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names : fourK, fiveK, sixK, sevenK, eightK, nineK, tenK, elevenK, twelveK, thirteenK, fourteenK, fifteenK, sixteenK, seventeenK, eighteenK, ...
min values : 5.611358e-07, 1.067142e-06, 2.508241e-06, 6.317553e-06, 2.286934e-05, 7.631922e-05, 1.272693e-04, 2.118215e-04, 2.602175e-04, 3.226203e-04, 4.390267e-04, 5.572032e-04, 7.313966e-04, 8.286005e-04, 8.297062e-04, ...
max values : 2.051069, 2.013452, 2.142908, 1.888403, 1.863014, 1.880628, 1.650615, 1.678033, 1.697732, 1.499115, 1.517264, 1.443677, 1.464867, 1.453581, 1.436394, ...
# Extract data to a matrix
Pop <- values(PopD.ALL)
r <- raster(PopD.ALL, 1)
r
class : RasterLayer
dimensions : 288, 720, 207360 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -60, 84 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : fourK
values : 5.611358e-07, 2.051069 (min, max)
We need to justify our decision to use a GAM over other models. This should include citations to back up those arguments.
We should make our decisions very transparent here. We should be able to justify our decision of 3 degrees of freedom over other possible values.
# Get the predctions from Population_trend script
load("prediction.RData")
# Read the polygons
origins <- readShapePoly('Origins_updated.shp')
# Extract data
cells <- do.call(rbind, sapply(per.origin, subset, select = 1))
#cells
g.means <- apply(prediction[-cells, ], 2, mean, na.rm = TRUE)
g.gams <- apply(prediction[-cells, ], 2, sd, na.rm = TRUE)
g.means2 <- apply(prediction[cells, ], 2, mean, na.rm = TRUE)
g.gams2 <- apply(prediction[cells, ], 2, sd, na.rm = TRUE)
#pdf("Global_pop_trend_comparisson.pdf", width = 25, height = 20)
par(mar = c(5, 7, 7, 5))
plot(seq(0, 1, length.out = length(time)) ~ time, col = "white", main = "GLOBAL",
xlim = c(21, 4), ylab = "Population Density (standardized)",
xlab = "Thousand of years ago", cex.lab = 1, cex.main = 1, cex.axis = 1)
down <- g.means - g.gams
up <- g.means + g.gams
lines(y = down, x = time, lty = 3, col = "gray40", lwd = 3)
lines(y = up, x = time, lty = 3, col = "gray40", lwd = 3)
lines(y = g.means, x = time, lwd = 4)
lines(y = g.means2, x = time, lwd = 3, col = "red")
down2 <- g.means2 - g.gams2
up2 <- g.means2 + g.gams2
lines(y = down2, x = time, lty = 3, col = "red", lwd = 3)
lines(y = up2, x = time, lty = 3, col = "red", lwd = 3)
polygon(cbind(c(12, 8.2, 8.2, 12, 12), c(-1, -1, 2, 2, -1)),
col = rgb(0, 1, 0, alpha = .2), border = F)
polygon(cbind(c(8.2, 4.2, 4.2, 8.2, 8.2), c(-1, -1, 2, 2, -1)),
col = rgb(.28, 0, .28, alpha = .2), border = F)
#dev.off()
# need to add a global mean, an everything but the origins mean, and a buffer around the origins mean.
# Read the polygons
origins <- readShapePoly('Origins_updated.shp')
# Extract data
per.origin <- extract(r, origins, cellnumber = TRUE, buffer = 100000)
names(per.origin) <- origins@data[, 1]
# Function standardization
std <- function(x) {
b <- (x - min(x)) / (max(x) - min(x))
return(rev(b))
}
# Calculating mean and
global.means <- global.SD <- list()
for (j in 1:length(per.origin)) {
#print(j)
originI <- Pop[per.origin[[j]][, 1], ]
time <- 21:4
originI <- na.exclude(originI)
b <- apply(originI, 1, std)
nJ <- nrow(originI)
predictions <- matrix(nrow = nJ, ncol = length(time))
colnames(predictions) <- as.character(time)
for(i in 1:nJ) {
# Need to show a gradient of these df values.
model <- gam(b[, i] ~ s(time, df = 15))
col <- sample(rainbow(100), 1)
predictions[i, ] <- predict(model)
}
global.means[[j]] <- apply(predictions, 2, mean)
global.SD[[j]] <- apply(predictions, 2, sd)
}
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
names(global.means) <- paste(names(per.origin), "Means")
names(global.SD) <- paste(names(per.origin), "SD")
plot(global.means[[1]], col=adjustcolor("cornflowerblue", alpha=0.8), pch=names(global.means[[1]]), type="b", xlab="year", ylab="Density", xaxt="n")
axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
global.means
$`W_African_Sav Means`
21 20 19 18 17 16 15
0.277236126 0.268146097 0.229965375 0.325239440 0.099371319 0.003780048 0.028853343
14 13 12 11 10 9 8
0.096533656 0.170551592 0.236171347 0.329881362 0.408474455 0.933326291 1.020506850
7 6 5 4
0.850021182 0.257169819 0.326479020 0.507805633
$`Sudanic_Savan Means`
21 20 19 18 17 16 15 14
0.01698931 0.01658303 0.01377887 0.01873437 0.02676287 0.05040199 0.11185374 0.17117563
13 12 11 10 9 8 7 6
0.22976620 0.29037257 0.39363331 0.46512887 0.95241742 1.01093274 0.89309987 0.35578816
5 4
0.54615181 0.75726423
$`West Africa T Means`
21 20 19 18 17 16 15 14
0.01898613 0.02119569 0.01844946 0.01229571 0.04011481 0.06425710 0.12281442 0.18024456
13 12 11 10 9 8 7 6
0.23876336 0.29694678 0.39395913 0.47594903 0.87575791 0.90858669 0.87451414 0.58149682
5 4
0.74507230 0.95664612
$`Ethipian plat Means`
21 20 19 18 17 16 15 14
0.09045898 0.09039249 0.08631005 0.08939636 0.10102942 0.11280965 0.16544499 0.23758729
13 12 11 10 9 8 7 6
0.33576370 0.42248406 0.56394916 0.62406635 0.94001810 0.95457968 0.92035401 0.56309708
5 4
0.60472053 0.91781126
$`NA Means`
21 20 19 18 17 16 15 14
0.66311422 0.66055581 0.64267837 0.68037338 0.57919737 0.44972219 0.29206295 0.19322360
13 12 11 10 9 8 7 6
0.13282795 0.11186128 0.07109712 0.13037398 0.67120744 0.78675575 0.73905383 0.83210726
5 4
0.65518007 0.51579662
$`E_North_Ameri Means`
21 20 19 18 17 16 15
0.008935466 0.010491576 0.013253014 0.005429808 0.033277010 0.067607455 0.151565641
14 13 12 11 10 9 8
0.252887930 0.364503617 0.473544250 0.693964954 0.767791047 0.692619560 0.694781027
7 6 5 4
0.698626071 0.419432398 0.858741718 1.009851858
$`New_Guinea Means`
21 20 19 18 17 16 15 14
0.28090841 0.27887032 0.26769422 0.28147012 0.25274529 0.21472830 0.16708455 0.11597030
13 12 11 10 9 8 7 6
0.08979978 0.08089207 0.01244318 0.04012938 0.56658193 0.60497997 0.64152040 0.95583067
5 4
0.70849108 0.66240742
$`Mesoamerica Means`
21 20 19 18 17 16 15
0.016195001 0.020624506 0.021438999 0.009346113 0.054411490 0.103081131 0.194265351
14 13 12 11 10 9 8
0.278924138 0.363943244 0.439908974 0.568842547 0.615137399 0.832150879 0.835544620
7 6 5 4
0.803694000 0.462771112 0.868405845 1.009287780
$`N_Lowland_SA Means`
21 20 19 18 17 16 15 14
0.19768919 0.19404181 0.19098494 0.19612871 0.19641766 0.18706802 0.16050661 0.13414511
13 12 11 10 9 8 7 6
0.10788427 0.09115094 0.02059283 0.04595208 0.59600575 0.62988681 0.65634027 0.95041179
5 4
0.71231938 0.76245058
$`NW_Lowland_SA Means`
21 20 19 18 17 16 15 14
0.30216048 0.30473388 0.30145828 0.30660886 0.29457158 0.26911208 0.24906510 0.20628929
13 12 11 10 9 8 7 6
0.15615286 0.12041864 0.03913096 0.03081728 0.29794031 0.39915852 0.50023797 0.94635445
5 4
0.69180450 0.60391276
$`Sava_W_India Means`
21 20 19 18 17 16 15
0.010624383 0.011415291 0.016124491 0.003918883 0.041128894 0.079911864 0.143474824
14 13 12 11 10 9 8
0.194472002 0.258964535 0.328023269 0.385019948 0.438200163 0.701906926 0.887863530
7 6 5 4
0.840257998 0.463636951 0.824692040 1.006026208
$`S_India Means`
21 20 19 18 17 16 15
0.011562958 0.018775299 0.021723311 0.007064695 0.040707283 0.071306979 0.139798009
14 13 12 11 10 9 8
0.225664613 0.281610839 0.348147871 0.403540117 0.477872218 0.800247791 0.852607259
7 6 5 4
0.838151683 0.653173598 0.915365428 1.005984490
$`Ganges_E_Indi Means`
21 20 19 18 17 16 15
0.011336515 0.012643677 0.013706701 0.007068004 0.034046963 0.061879159 0.107610274
14 13 12 11 10 9 8
0.154284669 0.214022369 0.261411689 0.354911485 0.405502151 0.701022088 0.856628658
7 6 5 4
0.857644757 0.634690322 0.846206514 1.003631222
$`Chinese_loess Means`
21 20 19 18 17 16 15
0.009851651 0.010480832 0.009006839 0.011611970 0.011877916 0.044124301 0.071626358
14 13 12 11 10 9 8
0.086101379 0.149418292 0.211519249 0.266988439 0.311227904 0.707572309 0.849429071
7 6 5 4
0.906202082 0.923045545 0.991676077 0.973623793
$`Japanese Means`
21 20 19 18 17 16 15 14
0.43296608 0.43672128 0.43035780 0.44830208 0.40355896 0.34383951 0.27578618 0.20037155
13 12 11 10 9 8 7 6
0.12621768 0.08676529 0.01311096 0.07301075 0.80925018 0.90329694 0.86470455 0.98192628
5 4
0.86398717 0.75827885
$`Lower-MiddleY Means`
21 20 19 18 17 16 15 14
0.01208663 0.01575036 0.01411414 0.01178427 0.03201081 0.05071188 0.09189183 0.12487305
13 12 11 10 9 8 7 6
0.15860810 0.20452035 0.24395821 0.28747052 0.68074944 0.77375375 0.84134458 0.69834173
5 4
0.90780965 1.00418953
$`South trop ch Means`
21 20 19 18 17 16 15 14
0.01671345 0.01747534 0.01513935 0.01802837 0.01568232 0.01227508 0.03863965 0.09161156
13 12 11 10 9 8 7 6
0.11464529 0.22358808 0.29758793 0.37545866 0.80524910 0.93801022 0.94552239 0.91083022
5 4
0.96600401 0.98801607
$`NA Means`
21 20 19 18 17 16 15 14
0.05066894 0.05230038 0.05292235 0.05835717 0.04889094 0.03936478 0.03905613 0.13386001
13 12 11 10 9 8 7 6
0.18460943 0.19718384 0.21980313 0.35557000 0.76574625 0.92168159 0.90805891 0.88692400
5 4
0.94897297 0.95173354
$`Southwes amaz Means`
21 20 19 18 17 16 15 14
0.19781479 0.19662496 0.18993796 0.19064183 0.20111408 0.19240561 0.16943207 0.14150350
13 12 11 10 9 8 7 6
0.11674148 0.10887194 0.06376160 0.07888967 0.60360527 0.61746807 0.65100738 0.92881641
5 4
0.92751249 0.92126285
$`C/S_Andes Means`
21 20 19 18 17 16 15 14
0.04013846 0.03976948 0.04037451 0.04009385 0.05985825 0.07147904 0.10113090 0.11532751
13 12 11 10 9 8 7 6
0.12894659 0.14465050 0.14431787 0.18643630 0.67781044 0.70463762 0.73107986 0.89570492
5 4
0.94294701 0.96109400
global.SD
$`W_African_Sav SD`
21 20 19 18 17 16 15 14
0.17269901 0.16622415 0.13867831 0.18544129 0.05180233 0.01767435 0.02327218 0.03336443
13 12 11 10 9 8 7 6
0.03541682 0.03519989 0.04853473 0.06084298 0.02424955 0.01252286 0.03085380 0.08163591
5 4
0.13119633 0.13670807
$`Sudanic_Savan SD`
21 20 19 18 17 16 15 14
0.01596981 0.01573248 0.01309909 0.01660033 0.01596607 0.02560278 0.04523896 0.05682794
13 12 11 10 9 8 7 6
0.06526664 0.06715609 0.05264207 0.05237042 0.01249360 0.02168305 0.03570080 0.14230242
5 4
0.12537535 0.13709822
$`West Africa T SD`
21 20 19 18 17 16 15 14
0.01673219 0.01578994 0.01481378 0.01136136 0.01969958 0.02756898 0.04194365 0.05049252
13 12 11 10 9 8 7 6
0.05967940 0.06586580 0.08288599 0.06990498 0.08183580 0.09246568 0.05351030 0.12807501
5 4
0.07954410 0.06756217
$`Ethipian plat SD`
21 20 19 18 17 16 15 14
0.16291219 0.16041160 0.15174222 0.17052526 0.12030900 0.08321087 0.10179650 0.15003746
13 12 11 10 9 8 7 6
0.20258904 0.22923829 0.26964589 0.25334530 0.04474472 0.05250856 0.04248866 0.18245586
5 4
0.07985206 0.10230497
$`NA SD`
21 20 19 18 17 16 15 14
0.31724393 0.31534358 0.30865208 0.32400309 0.29351537 0.25850263 0.20335450 0.15240433
13 12 11 10 9 8 7 6
0.09519069 0.06447755 0.08636095 0.12921181 0.14130818 0.18480720 0.14252473 0.14154360
5 4
0.09762597 0.11603513
$`E_North_Ameri SD`
21 20 19 18 17 16 15
0.008145198 0.008515398 0.009294010 0.004473189 0.011508708 0.015625265 0.022584915
14 13 12 11 10 9 8
0.025198449 0.024445301 0.022727993 0.032069814 0.040761082 0.052481734 0.052946139
7 6 5 4
0.039564858 0.037054913 0.026284034 0.001341543
$`New_Guinea SD`
21 20 19 18 17 16 15
0.082195833 0.082863879 0.078154136 0.082808803 0.078094798 0.071787935 0.053393900
14 13 12 11 10 9 8
0.049463057 0.038550008 0.027472188 0.011530095 0.008002951 0.087275661 0.089553441
7 6 5 4
0.056219416 0.005523272 0.033304801 0.044605179
$`Mesoamerica SD`
21 20 19 18 17 16 15
0.021153785 0.022939364 0.016641861 0.011819479 0.022523991 0.035286885 0.061362074
14 13 12 11 10 9 8
0.085539532 0.103613886 0.114220951 0.129846632 0.128105209 0.045201018 0.055611867
7 6 5 4
0.048259340 0.113311245 0.039613106 0.003286985
$`N_Lowland_SA SD`
21 20 19 18 17 16 15 14
0.10330655 0.10034704 0.10171417 0.10623824 0.09774895 0.08977896 0.06518589 0.04894999
13 12 11 10 9 8 7 6
0.03813627 0.02881662 0.02613161 0.02178818 0.11032995 0.12044717 0.10686664 0.01704901
5 4
0.08244668 0.12104669
$`NW_Lowland_SA SD`
21 20 19 18 17 16 15
0.130739697 0.132907733 0.132009694 0.133691638 0.127524837 0.113005201 0.069851646
14 13 12 11 10 9 8
0.057568162 0.059602563 0.063103075 0.072715681 0.066751796 0.078178507 0.093155849
7 6 5 4
0.058635826 0.005804702 0.045253754 0.052150366
$`Sava_W_India SD`
21 20 19 18 17 16 15
0.005786955 0.005220734 0.006138490 0.001834000 0.008370193 0.014809019 0.023478907
14 13 12 11 10 9 8
0.029202441 0.033611577 0.049744933 0.058186494 0.075529590 0.082953275 0.049622009
7 6 5 4
0.042997287 0.051914027 0.044540102 0.007729331
$`S_India SD`
21 20 19 18 17 16 15
0.011762488 0.015787105 0.011814528 0.006555100 0.012978839 0.012634400 0.017908110
14 13 12 11 10 9 8
0.048727505 0.042995471 0.028245014 0.035767644 0.056730798 0.038340653 0.034123821
7 6 5 4
0.031387785 0.026180677 0.027622224 0.001124058
$`Ganges_E_Indi SD`
21 20 19 18 17 16 15
0.009996104 0.010161509 0.010487678 0.006416099 0.014754648 0.023414633 0.037629153
14 13 12 11 10 9 8
0.051351876 0.064161762 0.070084828 0.099905287 0.108208542 0.147219682 0.125484765
7 6 5 4
0.102598416 0.092287701 0.093801548 0.002533703
$`Chinese_loess SD`
21 20 19 18 17 16 15
0.010457617 0.009801106 0.007833340 0.007565184 0.009443886 0.022651911 0.021857693
14 13 12 11 10 9 8
0.014895515 0.039283031 0.031720217 0.046021275 0.044369415 0.103392912 0.099384943
7 6 5 4
0.057457163 0.049823920 0.006320118 0.023000503
$`Japanese SD`
21 20 19 18 17 16 15
0.173874868 0.163342130 0.149019750 0.150323445 0.145901292 0.127512094 0.070015226
14 13 12 11 10 9 8
0.055198915 0.052044874 0.048228880 0.039058313 0.058440859 0.094423895 0.044348710
7 6 5 4
0.039955363 0.004581652 0.034395274 0.052817882
$`Lower-MiddleY SD`
21 20 19 18 17 16 15
0.012157572 0.013518960 0.012466976 0.011512224 0.014174312 0.016363557 0.021436642
14 13 12 11 10 9 8
0.032419006 0.039212090 0.043682275 0.056409379 0.056218411 0.050885376 0.046217941
7 6 5 4
0.029361739 0.048027193 0.040810241 0.002212724
$`South trop ch SD`
21 20 19 18 17 16 15 14
0.01251537 0.01380555 0.01211271 0.01426391 0.01218339 0.01093822 0.02713782 0.04214790
13 12 11 10 9 8 7 6
0.04932471 0.06819564 0.07256793 0.09153205 0.10740169 0.09418365 0.06698325 0.06470289
5 4
0.02305945 0.02213715
$`NA SD`
21 20 19 18 17 16 15 14
0.06537329 0.06279184 0.05844069 0.06515297 0.04831535 0.03065137 0.02933002 0.06506202
13 12 11 10 9 8 7 6
0.08780131 0.10016928 0.12958624 0.13599030 0.05745369 0.04854576 0.04360634 0.09354232
5 4
0.04283214 0.07353796
$`Southwes amaz SD`
21 20 19 18 17 16 15 14
0.10946970 0.11244845 0.11233937 0.11449979 0.09609846 0.07918164 0.04897765 0.04514773
13 12 11 10 9 8 7 6
0.05154662 0.06454150 0.08797545 0.08539320 0.04617412 0.04963616 0.04439196 0.07708328
5 4
0.04741189 0.06601495
$`C/S_Andes SD`
21 20 19 18 17 16 15 14
0.06693888 0.06209531 0.06430201 0.06473326 0.06066622 0.05122856 0.05624921 0.06267174
13 12 11 10 9 8 7 6
0.06421460 0.06696766 0.08404686 0.08570244 0.10382184 0.11311800 0.09010230 0.05474039
5 4
0.07111518 0.07379329
# Load patricks productivity PCA data
load('Productivity_ALL.RDATA')
# Load origin shapefiles
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
# Extract the data
prod.origin <- extract(Productivity.ALL, origins)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
# Mean and SD per region
means <- lapply(prod.origin, colMeans, na.rm = TRUE)
sds <- lapply(prod.origin, sd, na.rm = TRUE)
names(means) <- origins@data$CONTINENT
ymax <- max(unlist(means))
ymin <- min(unlist(means))
time <- 4:21
# Plot
#pdf("productivity.pdf", 20, 30)
par(mfrow = c(5, 4), mar = c(2, 2, 2, 0))
for (i in 1:length(means)) {
plot(y = means[[i]], x = time, xlim = c(21, 4), ylim = c(ymin, ymax),
main = names(means)[i], cex.main = 1, cex.lab = 1, cex.axis = 1,
ylab = "Productivity (PCA axis)", xlab = "Thousand of years ago (k)",
pch = 20, lwd = 1, type = "l",
col = c("purple", "green")[origin.time.region[i]])
up <- sds[[i]] + means[[i]]
down <- means[[i]] - sds[[i]]
lines(up ~ time, lty = 2)
lines(down ~ time, lty = 2)
}
#dev.off()
a <- layout(matrix(c(
1, 1, 1, 1, 1, 1, 1, 1,
3, 6, 7, 8, 9, 10, 11, 4,
3, 5, 5, 5, 5, 5, 5, 4,
3, 12, 13, 14, 15, 16, 17, 4,
2, 2, 2, 2, 2, 2, 2, 2
), 5, 8, byrow=TRUE), width=c(1, 1, 1, 1, 1, 1, 1, 1), height=c(0.5, 1, 1.5, 1, 0.5))
layout.show(a)
frameplot <- function(){
plot(21:0,rep(0, 22), xlim=c(17,4), ylim=c(0, 2.25), type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
frameplot_bottom <- function(){
plot(21:0,rep(0, 22), xlim=c(17,4), ylim=c(-0.25, 2), type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
frameplot()
frameplot_bottom()
d <- readPNG("earth.png")
png(file=paste("40962.png",sep=""),width=2000,height=1000, bg="transparent")
par(mar=c(0,0,0,0))
plot(seq(-180, 180, length.out = 19), seq(-90, 90, length.out = 19), type="n",xlim=c(-180, 180),ylim=c(-90, 90), xaxt="n")
rasterImage(d, -180, -90, 180, 90, interpolate=TRUE, col=d)
polygon(x=c(-180,-180, 180,180), y=c(-90, 90, 90, -90), col=adjustcolor("white", alpha=0.1))
#rasterImage(d, -13.5, -13.5, 375, 375, interpolate=TRUE, col=d)
plot(origins_subset, add=TRUE, col=adjustcolor("white", alpha=.8), xaxt="n", border="white", lwd=4) #still need to reproject!!!
dev.off()
###################
type_number <- 8
complex_figure <- function(type_number, i, means, sds){
if(i < 6) polygon(x=c(12,12,8.2,8.2), y=c(-1,3,3,-1), col=adjustcolor("cornflowerblue", alpha=0.4), border=NA)
if(i > 5) polygon(x=c(8.2,8.2,4.2,4.2), y=c(-1,3,3,-1), col=adjustcolor("limegreen", alpha=0.4), border=NA)
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(j)
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- -c(0, j(seq(0, maxer, length.out=100)))
#lines(x_seq, y_seq, type="l", ylim=c(-1,1))
#polygon(c(0, x_seq), c(0, y_seq), border="black", col=adjustcolor("cornflowerblue", alpha=0.5))
#abline(v= maxer - quantile(j)[2])
break_one_1 <- maxer
break_two_1 <- maxer - quantile(j)[2]
# polygon(x=c(break_two_1, break_two_1, break_one_1, break_one_1), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.5), border=NA)
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 10]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(j)
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 2+c(0, j(seq(0, maxer, length.out=100)))
#lines(x_seq, y_seq)
#polygon(c(0, x_seq), c(2, y_seq), border="black", col=adjustcolor("limegreen", alpha=0.5))
break_one_2 <- maxer
break_two_2 <- maxer - quantile(j)[2]
# polygon(x=c(break_two_2, break_two_2, break_one_2, break_one_2), y=c(1, 2, 2, 1), col=adjustcolor("limegreen", alpha=0.5), border=NA)
#abline(v=11)
type <- 1
if(type == 1){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
#polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
polygon(x=c(21:4,4:21), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
if(type == 2){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- x + 1 #scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
if(type == 3){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- x #scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
means_long_y <- c(1,1,1,1,1, meanss)
means_long_x <- c(0:4, 4:21)
break_one <- break_one_2
break_two <- break_two_2
# polygon(x=c(break_one, break_one, 22, 22), y=c(1, 2, 2, 1), col=adjustcolor("white", alpha=0.8), border=NA)
# polygon(x=c(break_two, break_two, break_one, break_one), y=c(1, 2, 2, 1), col=adjustcolor("white", alpha=0), border=NA)
# polygon(x=c(-1,-1, break_two , break_two), y=c(1.9, 3.1, 3.1, 1.9), col=adjustcolor("white", alpha=0.8), border=NA)
#abline(v= break_one, col="white")
#abline(v= break_two, col="white")
break_one <- break_one_1
break_two <- break_two_1
# polygon(x=c(break_one, break_one, 22, 22), y=c(0, 1, 1, 0), col=adjustcolor("white", alpha=0.8), border=NA)
# polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("white", alpha=0), border=NA)
# polygon(x=c(-1,-1, break_two , break_two), y=c(-1.1, .1, .1, -1.1), col=adjustcolor("white", alpha=0.8), border=NA)
#abline(v= break_one, col="white")
#abline(v= break_two, col="white")
#lines(x=c(break_one_2, break_one_2), y=c(1,3), col="white")
#lines(x=c(break_one_1, break_one_1), y=c(1,-1), col="white")
#lines(x=c(break_two_2, break_two_2), y=c(1,3), col="white")
#lines(x=c(break_two_1, break_two_1), y=c(1,-1), col="white")
#lines(4:21, meanss)
lines(21:4, meanss)
}
frameplot()
complex_figure(7, 1, global.means, global.SD)
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 3.5, 4, 4.5
axis(1)
axis(2)
quartz(width=8, height=8)
layout(matrix(c(
1, 1, 1, 1, 1, 1, 1, 1,
3, 6, 7, 8, 9, 10, 11, 4,
3, 5, 5, 5, 5, 5, 5, 4,
3, 12, 13, 14, 15, 16, 17, 4,
2, 2, 2, 2, 2, 2, 2, 2
), 5, 8, byrow=TRUE), width=c(1, 1, 1, 1, 1, 1, 1, 1), height=c(0.5, 1, 1.5, 1, 0.5))
par(mar=c(0,0,0,0))
# 1-4 label margins
blankplot <- function(){
plot(0,0, xlim=c(4,21), ylim=c(1, 1.25), bty="n", type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
blankplot()
blankplot()
blankplot()
blankplot()
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
as.character(origins$CONTINENT)
[1] "W_African_Sav" "Sudanic_Savan" "West Africa T" "Ethipian plat" "Fertile_Cresc"
[6] "E_North_Ameri" "New_Guinea" "Mesoamerica" "N_Lowland_SA" "NW_Lowland_SA"
[11] "Sava_W_India" "S_India" "Ganges_E_Indi" "Chinese_loess" "Japanese"
[16] "Lower-MiddleY" "South trop ch" "Chinese_loess" "Southwes amaz" "C/S_Andes"
#subset_order <- c(1, 2, 3, 5, 6, 8, 9, 10, 11, 12, 17, 18)
subset_order <- c(8, 10, 9, 5, 18, 7, 6, 20, 1, 2, 13, 16)
origins_subset <- origins[subset_order,]
origins_subset$CONTINENT
[1] Mesoamerica NW_Lowland_SA N_Lowland_SA Fertile_Cresc Chinese_loess New_Guinea
[7] E_North_Ameri C/S_Andes W_African_Sav Sudanic_Savan Ganges_E_Indi Lower-MiddleY
19 Levels: C/S_Andes Chinese_loess E_North_Ameri Ethipian plat ... West Africa T
d <- readPNG("earth.png")
png(file=paste("40962.png",sep=""),width=2000,height=1000, bg="transparent")
par(mar=c(0,0,0,0))
plot(seq(-180, 180, length.out = 19), seq(-90, 90, length.out = 19), type="n",xlim=c(-180, 180),ylim=c(-90, 90), xaxt="n")
rasterImage(d, -180, -90, 180, 90, interpolate=TRUE, col=d)
polygon(x=c(-180,-180, 180,180), y=c(-90, 90, 90, -90), col=adjustcolor("white", alpha=0.1))
#rasterImage(d, -13.5, -13.5, 375, 375, interpolate=TRUE, col=d)
plot(origins_subset, add=TRUE, col=adjustcolor("white", alpha=.8), xaxt="n", border="white", lwd=4) #still need to reproject!!!
dev.off()
quartz
2
d <- readPNG("40962.png")
dim(d)
[1] 1000 2000 4
par(mar=c(0,0,0,0))
plot(0:360,0:360,type="n",xlim=c(20,360),ylim=c(65,295), yaxt="n", xaxt="n")
rasterImage(d, -28.5, -13.5, 388, 375, interpolate=TRUE, col=d)
axis(2, label=seq(-90, 90, length.out = 19), at=seq(1, 360, length.out = 19), las=1)
mtext("latitude", 2, line=4, at=180)
abline(h=seq(1, 360, length.out = 19), col=adjustcolor("grey10", alpha= 0.4), lwd=1)
abline(h=180, col=adjustcolor("white", alpha= .5), lwd=1)
load('PopD_all_December.rdata')
# Extract the data
prod.origin <- extract(PopD.ALL, origins_subset)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
library(matrixStats)
# Mean and SD per region
means <- lapply(prod.origin, colMeans, na.rm = TRUE)
sds <- lapply(prod.origin, colSds, na.rm = TRUE)
## new values from Bruno's GAM model (produced in script called Population_Trend_per_y.R)
means <- global.means
sds <- global.gams
names(means) <- origins_subset@data$CONTINENT
ymax <- max(unlist(means))
ymin <- min(unlist(means))
time <- 4:21
#plot(origins)
#means[[1]] +
#sds[[1]]
#scale(as.numeric(means[[1]]), center=FALSE)
name_vector <- as.character(origins_subset@data$CONTINENT)
for(i in 1:12){
if(i > 6){frameplot()}else{frameplot_bottom()}
## customize polygons for each graph
if(i == 1){ #mesoamerica #values from Larson
complex_figure(3, i, means, sds)
}
#########
if(i == 2 ){ #NW lowlands SA #values from Larson
complex_figure(6, i, means, sds)
}
#########
if( i == 3){ #NW lowlands SA #values from Larson
complex_figure(6, i, means, sds)
}
#########
if(i == 4){ #Fertile crescent aka Southwest asia #values from Larson
complex_figure(8, i, means, sds)
}
#########
if(i == 5){ #loess plateau #values from Larson
complex_figure(2, i, means, sds)
}
#########
if(i == 6){ #new guinea #values from Larson
complex_figure(4, i, means, sds)
}
#########
if(i == 7){ #Eastern N.A. #values from Larson
complex_figure(5, i, means, sds)
}
#########
if(i == 8){ #Andes #values from Larson
complex_figure(6, i, means, sds)
}
#########
if(i == 9){ #W. African Sav #values from Larson
complex_figure(1, i, means, sds)
}
#########
if(i == 10){ #Sudanic sav #values from Larson
complex_figure(1, i, means, sds)
}
#########
if(i == 11){ #Ganges #values from Larson
complex_figure(7, i, means, sds)
}
#########
if(i == 12){ #loess #values from Larson
complex_figure(2, i, means, sds)
}
#lines(4:21, means[[i]])
abline(h = 1, col=adjustcolor("forestgreen", alpha=.5), lty=2)
# add axes to some locations
if(i == 1 | i == 7){axis(2, at=seq(0,2, by=0.25), label=seq(0,2, by=0.25), las=1)}
if(i == 6 | i == 12){axis(4, at=seq(0,2, by=0.25), label=seq(0,2, by=0.25), las=1)}
#if(i == 6 | i == 12){axis(4, at=seq(2,3, by=0.25), label=seq(0,1, by=0.25), las=1)
# axis(4, at=seq(-1,0, by=0.25), label=rev(seq(0,1, by=0.25)), las=1)
# }
if(i > 6){axis(1)} else{axis(3)}
# add text
if(i < 7){polygon(x=c(-30, -30, 30, 30), y=c(-0.1, -0.5, -0.5, -0.1), col="black")
mtext(name_vector[i], 1, line=-1.2, col="white", cex=0.5)}
if(i > 6){polygon(x=c(-30, -30, 30, 30), y=c(2.1, 2.5, 2.5, 2.1), col="black")
mtext(name_vector[i], 3, line=-1.2, col="white", cex=0.5)}
# add axis labels
if(i == 1 | i == 7){mtext("scaled density potential", 2, line=4, at=1)}
if(i == 3){mtext("Thousand years before present", 3, line=3.5, at =5)}
if(i == 9){mtext("Thousand years before present", 1, line=3.5, at =5)
}
}
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 4, 7, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:5] = 0, 1, 4, 7, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:6] = 0, 0.5, 1, ..., 2, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 0.5, 0.75, ..., 5, 7
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 0.5, 1, ..., 4.5, 7.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:1] = 0
Empirical CDF
Call: ecdf(maxer - match)
x[1:2] = 0, 4
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 0.5, 1.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 1, 1.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 2.25, 2.75, ..., 5.75, 6.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 2.25, 2.75, ..., 5.75, 6.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 3.5, 4, 4.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 0.5, 1, ..., 4.5, 7.5
saveToPDF <- function(...) {
d = dev.copy(pdf,...)
dev.off(d)
}
saveToPNG <- function(...) {
d = dev.copy(png,...)
dev.off(d)
}
## Try them out
saveToPDF("my.pdf", height=8,width=8)
quartz
2
saveToPNG("my.png", height=8, width=8, units="in", res=300)
quartz
2
dev.off()
null device
1